clc
clear
tic

% This script generates the tables and figures in the main text
%====================================================================%
%% =============== Figure 1: Illustrative Example ================= %%
%====================================================================%
gamma = 2;
nmu = 200;
mu_grid = linspace(0,1,nmu);
p = zeros(1,nmu);
for i = 1:nmu
    muA = mu_grid(i);
    p(i) = (1-muA)^(1/(gamma-1)); 
end
ones = ones(1,nmu);

figure
area(mu_grid,ones,'FaceColor',[0.7 0.7 0.7])
hold on
plot(mu_grid,p,'k','LineWidth',2)
area(mu_grid,p,'FaceColor',[0.3010 0.7450 0.9330])
xlabel('\mu_A')
ylabel('Matching Probability for Non-Outsourced Job (p)')
txt1 = 'Employment Effect Dominates';
text(0.15,0.25,txt1)
txt2 = 'Investment Disincentive Effect Dominates';
text(0.4,0.7,txt2)
yticks([0 0.2 0.4 0.6 0.8 1])
axis tight
saveas(gcf,'Figure1a','epsc')

clear
nmu = 200;
mu_grid = linspace(0,1,nmu);
gamma = 1.4;
p = zeros(1,nmu);
for i = 1:nmu
    muA = mu_grid(i);
    p(i) = (1-muA)^(1/(gamma-1)); 
end
ones = ones(1,nmu);

figure
area(mu_grid,ones,'FaceColor',[0.7 0.7 0.7])
hold on
plot(mu_grid,p,'k','LineWidth',2)
area(mu_grid,p,'FaceColor',[0.3010 0.7450 0.9330])
xlabel('\mu_A')
ylabel('Matching Probability for Non-Outsourced Job (p)')
txt1 = 'Employment Effect Dominates';
text(0.05,0.08,txt1)
txt2 = 'Investment Disincentive Effect Dominates';
text(0.4,0.5,txt2)
yticks([0 0.2 0.4 0.6 0.8 1])
axis tight
saveas(gcf,'Figure1b','epsc')

%====================================================================%
%% ======================= Specify Parameters ===================== %%
%====================================================================%
nu = 0.00622206;           % Quarterly probability of exiting the model
beta_discount = 0.99;      % Discount factor
beta = beta_discount*(1-nu); % Death-adjusted discount factor (factors in exit probability)

alpha =0.603158;           % Production function parameter
d =0.05316;                % Human capital depreciation rate
ell = 1.27;                % Matching function parameter
lambda =0.233887;          % Probability of search OTJ
chi = 3.782205;            % Investment cost parameter
b=0.692168;                % Unemployment benefit

muA =0.102465;             % Staffing agency fee
c =0.628833;               % Cost of posting a vacancy
delta =0.095524;           % Separation rate
h_sigma =0.200495;         % Standard deviation of initial human capital distribution
gamma=2;                   % Investment cost = (chi*i)^2

pr_phi(1) = 1.0-0.3084251; % Probability of entering the model with phi=1 (low fixed productivity type)
pr_phi(2) = 0.3084251;     % Probability of entering the model with phi=2 (high fixed productivity type)
z(1) = 1;                  % Productivity parameter f(phi,h) = z(phi)h^(alpha)
z(2) = 1.559533;

param = [nu beta alpha d ell lambda chi b muA c delta h_sigma pr_phi z];
param_orig = param;
nt = 400; 

% Specify grids %
nphi = 2;                  % Number of possible phi values (fixed productivity types)
hlb = 0.0;                 % Human capital grid lower bound 
hub = 7.0;                 % Human capital grid upper bound 
nh = 3000;                 % Number of possible human capital values
h_grid = linspace(hlb,hub,nh); % Create human capital grid

%--- Table 1: Assigned Parameters
param_assigned = [beta_discount ell nu pr_phi(2)];
param_calibrated = [lambda b muA c delta d chi alpha z(2) h_sigma];

%---- Data moments
[dm,pr_bachelors] = DO_combine_moments;          
rowNames = {'$\beta$','$\ell$','$\nu$','$\text{Pr}(\phi = H)$'};
description = {'Discount factor (quarterly)'
    'Matching function parameter'
    'Quarterly death probability'
    'Probability worker has Bachelor`s'};
values = {'0.99 (standard, 4\% annual interest)' '1.27 (den Haan et al., 2000)'
    '6.22206 $\times$ $10^{-3}$ (40 year career)' '0.30843 (ACS estimate)'};  

% Save as .tex file
fid = fopen('Table1.tex', 'w');
fprintf(fid, '\\begin{table}[H]\n\\centering\n');
fprintf(fid, '\\caption{Assigned Parameter Values}\n');
fprintf(fid, '\\begin{tabular}{c c c}\n');
fprintf(fid, '\\hline\n');
fprintf(fid, ' Parameter & Description & Value/Source \\\\\n');
fprintf(fid, '\\hline\n');
for i = 1:length(rowNames)
    fprintf(fid, '%s & %s & %s \\\\\n', rowNames{i}, description{i}, values{i});
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '\\label{tab: assigned_params}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%====================================================================%
%% ================== Solve Model and Report Moments ============== %%
%====================================================================%
indicator_policy=0; % Useful later when considering effects of policy changes
% Solve policy and value functions and steady-state distribution
[A_staffing,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,~,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          

% Plot implied total cost of creating each vacancy
cost_O = c+A_staffing;
cost_N = c;
figure
plot(h_grid,cost_O(1,:),'r','LineWidth',2)
hold on
plot(h_grid,cost_O(2,:),'-.b','LineWidth',2)
xlabel('Human Capital h (Relative to New Entrant Average=1)')
ylabel('Implied Total Cost of Creating Outsourced Match')
legend('Non-College','College','Location','NorthWest')
saveas(gcf,'Figure3a','epsc')

% Look at how expected benefit of entering either job type changes with c - across education types %
[t_Uhat_N_compare,t_Uhat_O_compare,c_grid,num_c] = DO_c_effect(param);

diff_O_minus_N = t_Uhat_O_compare - t_Uhat_N_compare;


% h=1
figure
plot(c_grid,reshape(diff_O_minus_N(1,430,:),[1,num_c]),'r','LineWidth',2)
hold on
plot(c_grid,reshape(diff_O_minus_N(2,430,:),[1,num_c]),'-.b','LineWidth',2)
yline(0,':k','LineWidth',1)
%xline(0.628833,'--k','LineWidth',1)
xlim([0.1 1])
ylabel('Outsourced Job Value Minus Exp. Not Outsourced Value')
xlabel('Firm Vacancy Cost c')
%title('h=1')
legend('Non-College (h=1)','College (h=1)','Location','NorthWest')
saveas(gcf,'Figure3b','epsc')

% Simulate wages 
simul_size = 20000;     % Number of individuals to follow in simulated sample
[simul,pr_entrant] = DO_sim_func(simul_size,param,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O);

% Generate Model Moments
[avg_diff_N,avg_diff_O,wage_sample_N1,wage_sample_N2,wage_sample_O1,wage_sample_O2,pr,pr_qrtr,pr_noncollege,pr_college,mm,wp_results,wp_noncollege_results,wp_college_results,wg_results,wg_noncollege,wg_college] = DO_moments_func(h_grid,simul,simul_size,param,dist_U,dist_N,dist_O,opt_ind_hp_N,opt_ind_hp_O,opt_p_U,opt_p_N,opt_p_O);
mm_orig = mm;
fprintf('Table 2: Jointly Calibrated Parameters - Data and Model Moments\n');
T = table(dm,mm,'VariableNames',{'Data','Model'},'RowName',{'Job-to-job transition rate','Unemployment Rate (%)','Percent Outsourced via staffing agency','Ratio outsourced with to w/out Bachelor`s Degree','Avg. Quarterly Separation Rate','Avg. Wage Loss After Unemployment Spell (%)','Wage Dispersion: p90/p50','Avg. Annual Real Wage Growth (%)',' College Wage Premium','Age 20 Unemployment (%)'});
disp(T)
mm_baseline=mm; % Save baseline model moments

%---- Table 2 ----%

rowNames = {'$\lambda$','$b_c$','$\mu^A$','$c$','$\delta$','$d$','$\chi$','$\alpha$', ...
            '$z(\phi_H)$','$\sigma_h$'};
fprintf('Table 2: Jointly Calibrated Parameters\n');
targets = ["Job-to-job transition rate";"Unemployment rate (\%)";"Percent outsourced via staffing agency";"Ratio outsourced w/out to outsourced w/ Bachelor's";"Average quarterly EU rate";"Avg. wage loss after unemployment spell (\%)";...
    "Wage dispersion: p90/p50";"Avg. annual real wage growth (\%)";"College wage premium";"Age 20 Unemployment rate (\%)"];
T = table(param_calibrated',targets,dm,mm,'VariableNames',{'Estimate','Targeted Moment','Data','Model'},'RowName',rowNames); 
disp(T)     

% Open file for writing LaTeX table
fid = fopen('Table2.tex', 'w');
fprintf(fid, '\\begin{table}[htbp]\n\\centering\n');
fprintf(fid, '\\caption{Jointly Calibrated Parameters}\n');
fprintf(fid, '\\scalebox{0.92}{%%\n');
fprintf(fid, '\\begin{tabular}{l c l c c}\n');
fprintf(fid, '\\hline\n');
fprintf(fid, 'Parameter & Estimate & Targeted Moment & Data & Model \\\\\n');
fprintf(fid, '\\hline\n');

% Loop through table rows
for i = 1:height(T)
    rowName  = T.Properties.RowNames{i};                 % Already contains LaTeX ($...$)
    estimate = T.Estimate(i);                            % numeric
    moment   = T.("Targeted Moment"){i};                 % string
    dataStr  = sprintf('%.3f', T.Data(i));               % format Data
    modelStr = sprintf('%.3f', T.Model(i));              % format Model
    fprintf(fid, '%s & %.6f & %s & %s & %s \\\\\n', ...
        rowName, estimate, moment, dataStr, modelStr);
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}}\n');
fprintf(fid, '\\label{tab: calibration}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%====================================================================%
%% ====================== Model Dynamics Results ================== %%
%====================================================================%
%-------- Wage Results
% Wage Penalty Results
wp_model(1,1) = wp_results(1);
wp_model(2,1) = wp_noncollege_results(1);
wp_model(3,1) = wp_college_results(1);

DO_ASEC_results = readtable('DO_ASEC_results.xlsx');

wp_data = table2array(DO_ASEC_results(1:3,4));  
wp_CIlb = table2array(DO_ASEC_results(1:3,5));
wp_CIub = table2array(DO_ASEC_results(1:3,6));

fprintf('Table 3: Wage Penalty Associated with Outsourcing\n');
T = table(wp_data,wp_model,'RowName',{'All Education Groups','No Bachelor`s Degree','Bachelor`s or More'},'VariableNames',{'Data Results','Model-Generated Data'});
disp(T)

% Row labels
row_labels = {'All Education Groups', 'No Bachelor`s Degree','Bachelor`s or More'};
% Open file
fileID = fopen('Table3.tex','w');
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Wage Penalty Associated with Outsourcing (Untargeted)}\n');
fprintf(fileID, '\\label{tab:untargeted_wp}\n');
fprintf(fileID, '\\scalebox{0.87}{\n');
fprintf(fileID, '\\begin{tabular}{ c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '  & Data Results & Model-Generated Data \\\\ \n');
fprintf(fileID, '\\hline\n');
for i = 1:3
    % Main data row
    fprintf(fileID, ' %s & %.3f & %.3f \\\\\n', row_labels{i}, wp_data(i), wp_model(i));
    % Confidence interval row
    fprintf(fileID, '       & [%.3f,%.3f] & \\\\\n', wp_CIlb(i), wp_CIub(i));
end
fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);

% Wage Growth Results
wg_model(1,1) = wg_results(1);
wg_model(2,1) = wg_noncollege(1);
wg_model(3,1) = wg_college(1);
wg_data = table2array(DO_ASEC_results(1:3,7));  
wg_CIlb = table2array(DO_ASEC_results(1:3,8));
wg_CIub = table2array(DO_ASEC_results(1:3,9));

fprintf('Table 4: Real Wage Growth Difference Outsourcing Coefficient\n');
T = table(wg_data,wg_model,'RowName',{'All Education Groups','No Bachelor`s Degree','Bachelor`s or More'},'VariableNames',{'Data Results','Model-Generated Data'});
disp(T)

row_labels = {'All Education Groups', 'No Bachelor`s Degree','Bachelor`s or More'};
fileID = fopen('Table4.tex','w');
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Annual Wage Growth Difference Associated with Outsourcing (Untargeted)}\n');
fprintf(fileID, '\\label{tab:untargeted_wg}\n');
fprintf(fileID, '\\scalebox{0.87}{\n');
fprintf(fileID, '\\begin{tabular}{ c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '  & Data Results & Model-Generated Data \\\\ \n');
fprintf(fileID, '\\hline\n');
for i = 1:3
    % Main data row
    fprintf(fileID, ' %s & %.3f & %.3f \\\\\n', row_labels{i}, wg_data(i), wg_model(i));
    % Confidence interval row
    fprintf(fileID, '       & [%.3f,%.3f] & \\\\\n', wg_CIlb(i), wg_CIub(i));
end
fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);

%---- Annual Transition Probabilities 
pr_data = table2array(DO_ASEC_results(:,1));
pr_data_noncollege = table2array(DO_ASEC_results(:,2));
pr_data_college = table2array(DO_ASEC_results(:,3));

pr = pr';
pr_noncollege = pr_noncollege';
pr_college = pr_college';

fprintf('Table 5: Annual Transition Probabilities (Untargeted)\n');
T = table(pr_data,pr,pr_data_noncollege,pr_noncollege,pr_data_college,pr_college,'VariableNames',{'Data_All_Groups','Model_All_Groups','Data_No_Bachelors','Model_No_Bachelors','Data_Bachelors','Model_Bachelors'},'RowName',{'Not_Outsourced_to_Unemployed','Outsourced_to_Unemployed','Not_Outsourced_to_Outsourced','Outsourced_to_Not_Outsourced','Unemployed_to_Not_Outsourced','Unemployed_to_Outsourced'});
disp(T)

% Write Table 5 to a LaTeX file
row_labels = {'Not Outsourced to Unemployed';'Outsourced to Unemployed';
    'Not Outsourced to Outsourced';'Outsourced to Not Outsourced';
    'Unemployed to Not Outsourced';'Unemployed to Outsourced'};
fileID = fopen('Table5.tex','w');

% Write LaTeX table preamble
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\begin{center}\n');
fprintf(fileID, '\\caption{Annual Transition Probabilities (Untargeted)}\n');
fprintf(fileID, '\\scalebox{0.87}{\n');
fprintf(fileID, '\\begin{tabular}{ c  c c  c c  c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '  & \\multicolumn{2}{c}{All Education Groups} & \\multicolumn{2}{c}{No Bachelor''s Degree}& \\multicolumn{2}{c}{Bachelor''s or More}\\\\\n');
fprintf(fileID, '      & Data & Model & Data & Model & Data & Model  \\\\ \n');
fprintf(fileID, '\\hline\n');
for i = 1:6
    fprintf(fileID, ' %s & %.4f & %.4f & %.4f & %.4f & %.4f & %.4f \\\\\n', ...
        row_labels{i}, ...
        pr_data(i), pr(i), ...
        pr_data_noncollege(i), pr_noncollege(i), ...
        pr_data_college(i), pr_college(i));
end
fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}}\n');
fprintf(fileID, '\\label{table:transition_probabilities}\\\\\n');
fprintf(fileID, '\\end{center}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);

%--- Quarterly Transition Probabilities Table (For Appendix) (Table 13)
pr_qrtr = pr_qrtr';
% Define row labels (transition descriptions)
row_labels = {'Not Outsourced to Unemployed','Outsourced to Unemployed','Not Outsourced to Outsourced', ...
    'Outsourced to Not Outsourced','Unemployed to Not Outsourced','Unemployed to Outsourced'};
fid = fopen('Table13.tex', 'w');
fprintf(fid, '\\begin{table}[h]\n');
fprintf(fid, '\\centering\n');
fprintf(fid, '\\caption{Annual and Quarterly Model Transition Probabilities}\n');
fprintf(fid, '\\begin{tabular}{ c c c }\n');
fprintf(fid, '\\hline\n');
fprintf(fid, '    & Annual & Quarterly \\\\\n');
fprintf(fid, '\\hline\n');
for i = 1:length(row_labels)
    fprintf(fid, ' %s & %.4f & %.4f \\\\\n', row_labels{i}, pr(i), pr_qrtr(i));
end
fprintf(fid, ' \\hline\n');
fprintf(fid, '\\end{tabular}\n');
fprintf(fileID, '\\label{tab:annual_qrtrly_transitions}\\\\\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);
%--- Entrant Transition Probabilities Table (Table 14)
fid = fopen('Table14.tex', 'w');
fprintf(fid, '\\begin{table}[h]\n');
fprintf(fid, '\\centering\n');
fprintf(fid, '\\caption{Entrant Quarterly Transitions}\n');
fprintf(fid, '\\begin{tabular}{ c c c }\n');
fprintf(fid, '\\hline\n');
fprintf(fid, ' Entrant to Unemployed & Entrant to Not Outsourced & Entrant to Outsourced \\\\\n');
fprintf(fid, '\\hline\n');
fprintf(fid, ' %.4f & %.4f & %.4f \\\\\n', pr_entrant(1), pr_entrant(2), pr_entrant(3));
fprintf(fid, '\n \\hline\n');
fprintf(fid, '\\end{tabular}\n');
fprintf(fileID, '\\label{tab:entrant_transitions}\\\\\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%-------------- Model Selection into Outsourcing --------------------%
perc_college_N = (( sum(dist_N(2,:)))/( sum(sum(dist_N(:,:))) ))*100;
perc_college_O = (( sum(dist_O(2,:)))/( sum(sum(dist_O(:,:))) ))*100;

% Selection: Average h
avg_h_N = 0;
avg_h_O = 0;
for ind_h = 1:nh
    h = h_grid(ind_h);

    avg_h_N = avg_h_N + h*(sum(dist_N(:,ind_h)) );
    avg_h_O = avg_h_O + h*(sum(dist_O(:,ind_h)) );
end
avg_h_N = avg_h_N/( sum(sum(dist_N(:,:))) );
avg_h_O = avg_h_O/( sum(sum(dist_O(:,:))) );

fprintf('Table 6: Model Selection Into Outsourcing\n');
T = table([perc_college_O;avg_h_O],[perc_college_N;avg_h_N],'VariableNames',{'Outsourced','Non_Outsourced'},'RowName',{'Percent_with_Bachelors','Average_h'});
disp(T)

% Write Table 6 to LaTeX file
fileID = fopen('Table6.tex','w');
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Model Selection Into Outsourcing}\n');
fprintf(fileID, '\\label{tab:Selection}\n');
fprintf(fileID, '\\scalebox{0.9}{\n');
fprintf(fileID, '\\begin{tabular}{ c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '  & Outsourced & Non-Outsourced \\\\ \n');
fprintf(fileID, '\\hline\n');

fprintf(fileID, ' Percent with Bachelor''s ($\\phi=H$) & %.3f\\%% & %.3f\\%% \\\\\n', perc_college_O, perc_college_N);
fprintf(fileID, ' Average h (Relative to New Entrant Average = 1) & %.3f & %.3f \\\\\n', avg_h_O, avg_h_N);

fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);
%------------------ Figure 4: Optimal Search Decisions ------------------------%
figure
plot(h_grid,opt_p_U(1,:),'r','LineWidth',2)
hold on
plot(h_grid,opt_p_U(2,:),'-.b','LineWidth',2)
xlim([0.0 5])
legend('Non-College','College','Location','NorthEast')
xlabel('Human Capital h (Relative to New Entrant Average=1)')
ylabel('Optimal Worker Matching Probability (p)')
saveas(gcf,'Figure4a','epsc')

figure
plot(h_grid,opt_p_O(1,:),'r','LineWidth',2)
hold on
plot(h_grid,opt_p_O(2,:),'-.b','LineWidth',2)
xlim([0.0 5])
legend('Non-College','College','Location','NorthWest')
xlabel('Human Capital h (Relative to New Entrant Average=1)')
ylabel('Optimal Worker Matching Probability (p)')
saveas(gcf,'Figure4b','epsc')

%--------------- Figure 5: Human Capital Investment and Distributions -----------------%
opt_hp_N = h_grid(opt_ind_hp_N);
opt_hp_O = h_grid(opt_ind_hp_O);

figure
plot(h_grid,opt_hp_N(1,:),':','LineWidth',2)
hold on
plot(h_grid,opt_hp_N(2,:),'Color',[0.9290 0.6940 0.1250],'LineWidth',2)
plot(h_grid,h_grid,':k','LineWidth',1)
xlim([0.0 5])
yticks([0 1 2 3 4 5])
legend('Non-College','College','45-Degree Line','Location','NorthWest')
xlabel('Current Human Capital h (Relative to New Entrant Mean=1)')
ylabel('Optimal Next Period Human Capital (Non-Outsourced Job)')
saveas(gcf,'Figure5a','epsc')

% HC distributions
cdf_N= zeros(2,nh);
cdf_O= zeros(2,nh);
for ind_h =1:nh
    for ind_phi=1:nphi
        if (ind_h==1)
            cdf_N(ind_phi,ind_h) = (dist_N(ind_phi,ind_h))/( sum(dist_N(ind_phi,:)));
            cdf_O(ind_phi,ind_h) = (dist_O(ind_phi,ind_h))/( sum(dist_O(ind_phi,:)));
        else
            cdf_N(ind_phi,ind_h) = cdf_N(ind_phi,ind_h-1) + (dist_N(ind_phi,ind_h))/( sum(dist_N(ind_phi,:)) );
            cdf_O(ind_phi,ind_h) = cdf_O(ind_phi,ind_h-1) + (dist_O(ind_phi,ind_h))/( sum(dist_O(ind_phi,:)) );
        end
    end
end

figure
plot(h_grid,cdf_N(1,:),'LineStyle',':','Color',[0 0.4470 0.7410],'LineWidth',2)
hold on
plot(h_grid,cdf_O(1,:),'Color',[0.8500 0.3250 0.0980],'LineWidth',2)
plot(h_grid,cdf_N(2,:),'LineStyle','--','Color',[0.9290 0.6940 0.1250],'LineWidth',2)
plot(h_grid,cdf_O(2,:),'LineStyle','-.','Color',[0.4940 0.1840 0.5560],'LineWidth',2)
legend('Non-College: Not Outsourced','Non-College: Outsourced','College: Not Outsourced','College: Outsourced','location','NorthWest')
ylabel('CDF')
xlabel('Human Capital (Relative to New Entrant Mean Normalized to 1)')
ylim([0 1])
xlim([0.0 5])
saveas(gcf,'Figure5b','epsc')

%---------------- Figure 6: Optimal Investment Rule: h'-h ----------------%
h_groups = linspace(1,length(avg_diff_N(1,:)),length(avg_diff_N(1,:)));
figure
scatter(h_groups,avg_diff_N(1,:),'b','filled')
hold on
plot(h_groups,avg_diff_N(1,:),':b','HandleVisibility','off')
scatter(h_groups,avg_diff_O(1,:),'s','r','filled')
plot(h_groups,avg_diff_O(1,:),':r','HandleVisibility','off')
scatter(h_groups,avg_diff_N(2,:),'>','k','filled')
plot(h_groups,avg_diff_N(2,:),':k','HandleVisibility','off')
scatter(h_groups,avg_diff_O(2,:),'x','s','MarkerEdgeColor',[0.9290 0.6940 0.1250],'MarkerFaceColor',[0.9290 0.6940 0.1250])
plot(h_groups,avg_diff_O(2,:),'Color',[0.9290 0.6940 0.1250],'LineStyle',':','HandleVisibility','off')
yline(0,'k')
legend('Non-College: Not Outsourced','Non-College: Outsourced','College: Not Outsourced','College: Outsourced')
ylabel('Optimal Human Capital Increase: h`-h')
xlabel('Current Human Capital (h)')
xticks([1 3 5 7 9 11 ])
xticklabels({'[0,0.25)','[0.5,0.75)','[1,1.25)','[1.5,1.75)','[2,2.25)','[2.5,2.75)'})
axis tight
box on
saveas(gcf,'Figure6','epsc')

%====================================================================%
%% ======== 5.3 Counterfactual Economy without Outsourcing ======== %%
%====================================================================%

% Loop over "policy" 
% policy = 1 means the original economy with outsourcing - save moments of interest %
% policy = 2 means the economy without outsourcing - compute new steady state and save moments of interest %

npolicy = 2; % Number of policies considered
results = zeros(18,npolicy);
avg_hc_sim = zeros(npolicy,nt);

for policy = 1:npolicy
    
    if policy==2 % Need to re-solve the model in this case
        muA=1;
        param = [nu beta alpha d ell lambda chi b muA c delta h_sigma pr_phi z];
        % Solve policy and value functions and steady-state distribution
        [A_staffing,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,~,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          

        % Simulate wages 
        [simul,~] = DO_sim_func(simul_size,param,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O);
    end
    
    % Compute and save moments of interest for comparison
    %-- Unemployment rate
    results(1,policy) = (sum(sum(dist_U(:,:)))/( sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:))) + sum(sum(dist_U(:,:))) ))*100;
    results(2,policy) = (sum(dist_U(2,:))/( sum(dist_N(2,:)) + sum(dist_O(2,:)) + sum(dist_U(2,:)) ))*100;
    results(3,policy) = (sum(dist_U(1,:))/( sum(dist_N(1,:)) + sum(dist_O(1,:)) + sum(dist_U(1,:)) ))*100;

    %----------- Output from Employment (GDP)
    for ind_phi = 1:nphi
        for ind_h = 1:nh
            h = h_grid(ind_h);
            results(4,policy) = results(4,policy) + (z(ind_phi)*h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
            if (ind_phi==2)
                results(5,policy) = results(5,policy) + (z(ind_phi)*h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
            else
                results(6,policy) = results(6,policy) + (z(ind_phi)*h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
            end
        end
    end

    %---------- Average Human Capital
    h_sum_total = 0;
    h_sum_nc = 0;
    h_sum_coll = 0;

    for ind_phi = 1:nphi
        for ind_h = 1:nh
            h = h_grid(ind_h);

            results(7,policy) = results(7,policy) + h*(dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
            h_sum_total = h_sum_total + dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h);

            if (ind_phi==2)
                results(8,policy) = results(8,policy) + h*(dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
                h_sum_coll = h_sum_coll + dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h);
            else
                results(9,policy) = results(9,policy) + h*(dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
                h_sum_nc = h_sum_nc + dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h);
            end
        end
    end
    results(7,policy) = results(7,policy)/h_sum_total;
    results(8,policy) = results(8,policy)/h_sum_coll;
    results(9,policy) = results(9,policy)/h_sum_nc;
   
    %-------- Average Wage
    count_wage_total = 0;
    count_wage_coll = 0;
    count_wage_nc = 0;

    for id = 1:simul_size
        for t = 1:nt
            if (simul(id,t).stat~=999 && simul(id,t).stat~=0) % If employed
                results(10,policy)  = results(10,policy) + simul(id,t).wage;
                count_wage_total = count_wage_total + 1;
                if (simul(id,t).phi_ind==2)
                    results(11,policy) = results(11,policy) + simul(id,t).wage;
                    count_wage_coll = count_wage_coll + 1;
                else
                    results(12,policy) = results(12,policy) + simul(id,t).wage;
                    count_wage_nc = count_wage_nc + 1;
                end
            end
        end
    end
    results(10,policy) = results(10,policy)/count_wage_total;
    results(11,policy) = results(11,policy)/count_wage_coll;
    results(12,policy) = results(12,policy)/count_wage_nc;
    
    %--------- Total Utility
    util = 0;
    for ind_phi = 1:nphi
        for ind_h = 1:nh
            h = h_grid(ind_h);
            hp_N = h_grid(opt_ind_hp_N(ind_phi,ind_h));
            hp_O = h_grid(opt_ind_hp_O(ind_phi,ind_h));
            % Compute production minus investment cost
            util = util + dist_U(ind_phi,ind_h)*(b*z(ind_phi)*(h^alpha))...
                + z(ind_phi)*(h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h))...
                - (((hp_N - h*(1-d))*chi)^gamma)*dist_N(ind_phi,ind_h) - (((hp_O - h*(1-d))*chi)^gamma)*dist_O(ind_phi,ind_h);
            % Now subtract out vacancy posting costs
            if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3) % Search in unemployment for non-outsourced job
                if (~isnan(opt_theta_U(ind_phi,ind_h)))
                    util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*c;
                end
            else                                              % Search in unemployment for outsourced job
                if (~isnan(opt_theta_U(ind_phi,ind_h)))
                    util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*(c + A_staffing(ind_phi,ind_h));
                end
            end

            if (GN(ind_phi,ind_h)==1 || GN(ind_phi,ind_h)==3) % Search from non-outsourced job for non-outsourced job
                if (~isnan(opt_theta_N(ind_phi,ind_h)))
                    util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*c;
                end
            else                                              % Search from outsourced job for non-outsourced job
                if (~isnan(opt_theta_N(ind_phi,ind_h)))
                    util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                end
            end

            if (GO(ind_phi,ind_h)==1 || GO(ind_phi,ind_h)==3) % Search from outsourced job for non-outsourced job
                if (~isnan(opt_theta_O(ind_phi,ind_h)))
                    util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*c;
                end
            else
                if (~isnan(opt_theta_O(ind_phi,ind_h)))
                    util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                end
            end
        end
    end
    results(13,policy) = util;
    
    % Now compute consumption only among non-college or college
    for ind_phi= 1:nphi
        util = 0;
        for ind_h = 1:nh
            h = h_grid(ind_h);
            hp_N = h_grid(opt_ind_hp_N(ind_phi,ind_h));
            hp_O = h_grid(opt_ind_hp_O(ind_phi,ind_h));

            util = util + dist_U(ind_phi,ind_h)*(b*z(ind_phi)*(h^alpha))...
                + z(ind_phi)*(h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h))...
        	    - (((hp_N - h*(1-d))*chi)^gamma)*dist_N(ind_phi,ind_h) - (((hp_O - h*(1-d))*chi)^gamma)*dist_O(ind_phi,ind_h);
            
            if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)
                if (~isnan(opt_theta_U(ind_phi,ind_h)))
                    util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*c;
                end
            else
                if (~isnan(opt_theta_U(ind_phi,ind_h)))
                    util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*(c + A_staffing(ind_phi,ind_h));
                end
            end

            if (GN(ind_phi,ind_h)==1 || GN(ind_phi,ind_h)==3)
                if (~isnan(opt_theta_N(ind_phi,ind_h)))
                    util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*c;
                end
            else
                if (~isnan(opt_theta_N(ind_phi,ind_h)))
                    util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                end
            end

            if (GO(ind_phi,ind_h)==1 || GO(ind_phi,ind_h)==3)
                if (~isnan(opt_theta_O(ind_phi,ind_h)))
                    util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*c;
                end
            else
                if (~isnan(opt_theta_O(ind_phi,ind_h)))
                    util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                end
            end
        end
        if (ind_phi==2)
            results(14,policy) = util;  
        else
            results(15,policy) = util; 
        end
    end % end loop over phi possibilities

    %-------- Entrant Value
    %---------- New Entrant Value
    results(16,policy) = pr_phi(1)*U_hat(1,1) + pr_phi(2)*U_hat(2,1);
    results(17,policy) = U_hat(2,1);
    results(18,policy) = U_hat(1,1);

    % With simulated data "simul" compare average human capital by time spent in the model t %
    for t = 1:nt
        count_avgh_t = 0;
        for id = 1:simul_size
            if (simul(id,t).stat~=999)
                count_avgh_t = count_avgh_t + 1;
                avg_hc_sim(policy,t) = avg_hc_sim(policy,t) + h_grid(simul(id,t).h_ind);
            end
        end
        avg_hc_sim(policy,t) = avg_hc_sim(policy,t)/count_avgh_t;
    end

end % end loop over policy

disp_results = zeros(18);
for j=1:18
    if (j<4) % These are the unemployment rate results, where percentage point change is reported
        disp_results(j) = results(j,2)-results(j,1);
    else     % These are all non-unemployment rate results, where percent change is reported
        disp_results(j) = ((results(j,2)-results(j,1))/results(j,1))*100;
    end
end

total_results = [disp_results(1); disp_results(4); disp_results(7); disp_results(10); disp_results(13); disp_results(16)];
coll_results = [disp_results(2); disp_results(5); disp_results(8); disp_results(11); disp_results(14); disp_results(17)];
nc_results = [disp_results(3); disp_results(6); disp_results(9); disp_results(12); disp_results(15); disp_results(18)]; 

%------- Table 7 ------%
fprintf('Table 7: Steady State Comparison: General Equilibrium Effect of Removing Outsourcing\n');
T = table(total_results,nc_results,coll_results,'RowName',{'Unemployment','Output','Average_HC','Average_Wage','Total_Utility','Entrant_Value'},'VariableNames',{'Full_Population','Only_Non_College','Only_College'});
disp(T)

% Write table 7 to a LaTeX file
row_labels = {'Unemployment rate';'Output from employment';'Average human capital';
    'Average wage';'Total consumption';'Entrant value'};

% Units for each row (e.g., "pp" or "\%%" for LaTeX)
row_units = {'pp';'\%';'\%';'\%';'\%';'\%'};

fileID = fopen('Table7.tex','w');
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Steady State Comparison: Aggregate Effects of Removing Staffing Agencies}\n');
fprintf(fileID, '\\label{tab:Counterfactual_No_Outsourcing}\n');
fprintf(fileID, '\\begin{tabular}{ c c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '     & Full Population & Only Non-College & Only College \\\\ \n');
fprintf(fileID, '\\hline\n');

for i = 1:6
    % Print with sign (+/-), 3 decimal places, and units
    fprintf(fileID, ' %s & %+.3f%s & %+.3f%s & %+.3f%s \\\\\n', ...
        row_labels{i}, ...
        total_results(i), row_units{i}, ...
        nc_results(i), row_units{i}, ...
        coll_results(i), row_units{i});
end
fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);


% With simulated data "simul" compare average human capital by time spent in the model t %
% First compute % change in avg hc after the policy change 
pc_avg_hc = zeros (1,nt);
for t=1:nt
    pc_avg_hc(t) = ((avg_hc_sim(2,t) - avg_hc_sim(1,t))/avg_hc_sim(1,t))*100;
end

t_grid = linspace(1,nt,nt);
figure
plot(t_grid,pc_avg_hc,'r','LineWidth',2)
hold on
yline(0,':k')
ylabel('Percent Change in Simulated Average Human Capital')
yticks([-6 -5 -4 -3 -2 -1 0 1])
yticklabels({'-6%','-5%','-4%','-3%','-2%','-1%','0%','1%'})
xticks([1 41 81 121 161 201])
xticklabels({'20','30','40','50','60','70'})
xlim([1 221])
xlabel('Age')
box on
saveas(gcf,'Figure7','epsc')
%====================================================================%
%% ======= 5.4 Employment and Investment Disincentive Channels in Full Quantitative Model ======= %%
%====================================================================%
% ------- For current c, loop over different muA values to create Table 12
muA_vals = linspace(0.1,0.4,7);
b_vals = [0.692168 0.02465	-0.7173	-1.539	-2.459	-3.479	-4.602];
results = zeros(2,7);
urate_54 = zeros(2,7);
for ind_params = 1:7

    npolicy = 2; % Number of policies considered - Original steady state and economy without outsourcing %
    avg_h_total = zeros(1,npolicy);
    util_total = zeros(1,policy);
    for policy = 1:2 

        if policy==1
            param = param_orig;
            if (ind_params==1)
            else
                param(9) = muA_vals(ind_params);  % muA
            end
            param(8) = b_vals(ind_params);        % b

            indicator_policy=0;
            % Solve policy and value functions and steady-state distribution
            [A_staffing,~,~,~,~,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,~,~,~,~,~,~,~,~,~,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          
        else
            % Need to re-solve the model in this case
            param = param_orig;
            param(9) =1;                           % muA 
            param(8) = b_vals(ind_params);         % b               
            indicator_policy=0;
            % Solve policy and value functions and steady-state distribution
            [A_staffing,~,~,~,~,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,~,~,~,~,~,~,~,~,~,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          
        end
    
        % Compute and save moments of interest for comparison
        %---------- Average Human Capital
        h_sum_total = 0;
        for ind_phi = 1:nphi
            for ind_h = 1:nh
                h = h_grid(ind_h);

                avg_h_total(policy)  = avg_h_total(policy)  + h*(dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
                h_sum_total = h_sum_total + dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h);
            end
        end
        avg_h_total(policy) = avg_h_total(policy)/h_sum_total;

        %--------- Unemployment rate
        urate_54(policy,ind_params) = (sum(sum(dist_U(:,:)))/( sum(sum(dist_U(:,:))) + sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:)))) )*100;

        %--------- Total Utility
        util = 0;
        for ind_phi = 1:nphi
            for ind_h = 1:nh
                h = h_grid(ind_h);
                hp_N = h_grid(opt_ind_hp_N(ind_phi,ind_h));
                hp_O = h_grid(opt_ind_hp_O(ind_phi,ind_h));
                % Compute production minus investment cost
                util = util + dist_U(ind_phi,ind_h)*(b*z(ind_phi)*(h^alpha))...
                    + z(ind_phi)*(h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h))...
                    - (((hp_N - h*(1-d))*chi)^gamma)*dist_N(ind_phi,ind_h) - (((hp_O - h*(1-d))*chi)^gamma)*dist_O(ind_phi,ind_h);
                % Now subtract out vacancy posting costs
                if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3) % Search in unemployment for non-outsourced job
                    if (~isnan(opt_theta_U(ind_phi,ind_h)))
                        util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*c;
                    end
                else                                              % Search in unemployment for outsourced job
                    if (~isnan(opt_theta_U(ind_phi,ind_h)))
                        util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*(c + A_staffing(ind_phi,ind_h));
                    end
                end

                if (GN(ind_phi,ind_h)==1 || GN(ind_phi,ind_h)==3) % Search from non-outsourced job for non-outsourced job
                    if (~isnan(opt_theta_N(ind_phi,ind_h)))
                        util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*c;
                    end
                else                                              % Search from outsourced job for non-outsourced job
                    if (~isnan(opt_theta_N(ind_phi,ind_h)))
                        util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                    end
                end

                if (GO(ind_phi,ind_h)==1 || GO(ind_phi,ind_h)==3) % Search from outsourced job for non-outsourced job
                    if (~isnan(opt_theta_O(ind_phi,ind_h)))
                        util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*c;
                    end
                else
                    if (~isnan(opt_theta_O(ind_phi,ind_h)))
                        util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                    end
                end
            end
        end
        util_total(policy) = util;
    end % end loop over policy
    results(1,ind_params) = ((avg_h_total(2) - avg_h_total(1))/avg_h_total(1))*100;
    results(2,ind_params) = ((util_total(2) - util_total(1))/util_total(1))*100;
end % end loop over parameter values

results_urate_change = urate_54(2,:)-urate_54(1,:);


%-------- Table 8 

fileID = fopen('Table8.tex','w');

fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Effects of Removing Staffing Agencies}\n');
fprintf(fileID, '\\label{tab:Channels_Full_Model}\n');
fprintf(fileID, '\\scalebox{0.8}{\n');
fprintf(fileID, '\\begin{tabular}{ c c c c c c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, ' $\\mu_A$ & 0.102 (Table 2) & 0.15 & 0.20 & 0.25 & 0.30 & 0.35 & 0.40 \\\\ \n');
fprintf(fileID, '  $c$ & 0.629 (Table 2)  & 0.55 & 0.50 & 0.45 & 0.40 & 0.35 & 0.30 \\\\\n');
fprintf(fileID, '\\hline\n');

fprintf(fileID, ' Unemployment rate &');
for j = 1:7
    val = results_urate_change(j);
    if val >= 0
        fprintf(fileID, ' +%.3fpp', val);
    else
        fprintf(fileID, ' %.3fpp', val);
    end
    if j < 7
        fprintf(fileID, ' &');
    else
        fprintf(fileID, ' \\\\\n');
    end
end
fprintf(fileID, ' Avg. human capital &');
for j = 1:7
    val = results(1,j);
    if val >= 0
        fprintf(fileID, ' +%.3f\\%%', val);
    else
        fprintf(fileID, ' %.3f\\%%', val);
    end
    if j < 7
        fprintf(fileID, ' &');
    else
        fprintf(fileID, ' \\\\\n');
    end
end
fprintf(fileID, ' Total Utility &');
for j = 1:7
    val = results(2,j);
    if val >= 0
        fprintf(fileID, ' +%.3f\\%%', val);
    else
        fprintf(fileID, ' %.3f\\%%', val);
    end
    if j < 7
        fprintf(fileID, ' &');
    else
        fprintf(fileID, ' \\\\\n');
    end
end

fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}}\n\n');

% Add tablenotes section
fprintf(fileID, '\\begin{tablenotes}\n');
fprintf(fileID, '\\item \\footnotesize{\\textit{Notes:} When altering $\\mu^A$ and $c$, solve for the value of $b_c$ that keeps the percentage of workers in outsourced jobs constant in the baseline economy before the policy change.}\n');
fprintf(fileID, '\\end{tablenotes}\n\n');

fprintf(fileID, '\\end{table}\n');
fclose(fileID);



% ------- Loop over different muA and c values to create figure 8
muA_vals = linspace(0.1,0.4,7);
c_vals = linspace(0.6,0.3,7); 
% Input b values that keep % outourced constant - see Table 8
b_vals = [0.7106	0	-0.76	-1.605	-2.556	-3.611	-4.771;0.685	-0.053	-0.844	-1.738	-2.746	-3.869	-5.11;0.6546	-0.115	-0.948	-1.897	-2.974	-4.185	-5.52;0.618	-0.191	-1.078	-2.095	-3.263	-4.57	-6.027;0.571	-0.288	-1.235	-2.344	-3.618	-5.057	-6.663;0.5108	-0.408	-1.443	-2.665	-4.078	-5.687	-7.486;0.429	-0.567	-1.721	-3.098	-4.7	-6.538	-8.6];

hc_change = zeros(7,7);

for ind_c = 1:7
    for ind_muA = 1:7
  
        % Loop over "policy" 
        % policy = 1 means the original economy with outsourcing - save moments of interest %
        % policy = 2 means the economy without outsourcing - compute new steady state and save moments of interest %
        npolicy = 2; % Number of policies considered - Original steady state and economy without outsourcing %
        avg_h_total = zeros(1,npolicy);
        for policy = 1:2 
    
            if policy==1
                param = param_orig;
                param(9) = muA_vals(ind_muA);         % muA
                param(8) = b_vals(ind_c,ind_muA);     % b
                param(10) = c_vals(ind_c);            % c  

                indicator_policy=0;
                % Solve policy and value functions and steady-state distribution
                [~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          
            else
                % Need to re-solve the model in this case
                param = param_orig;
   
                param(8) = b_vals(ind_c,ind_muA);       % b
                param(10) = c_vals(ind_c);              % c  

                param(9) =1;                            % muA 
                indicator_policy=0;
                % Solve policy and value functions and steady-state distribution
                [~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          

            end
    
            % Compute and save moments of interest for comparison
    
            %---------- Average Human Capital
            h_sum_total = 0;
            for ind_phi = 1:nphi
                for ind_h = 1:nh
                    h = h_grid(ind_h);

                    avg_h_total(policy)  = avg_h_total(policy)  + h*(dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
                    h_sum_total = h_sum_total + dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h);
                end
            end
            avg_h_total(policy) = avg_h_total(policy)/h_sum_total;

        end % end loop over policy
        hc_change(ind_c,ind_muA) = ((avg_h_total(2) - avg_h_total(1))/avg_h_total(1))*100;

    end % end loop over muA
end % end loop over c

%==== Now Plot Results ======%
muA_grid = linspace(0, 6, 7); 
c_grid = linspace(0, 6, 7)';   
% Create meshgrid
[muA_mesh, c_mesh] = meshgrid(muA_grid, c_grid);

% Plot the surface
figure
s = surf(muA_mesh, c_mesh, hc_change);
shading interp
colormap('parula')
%colorbar
s.FaceAlpha = 0.7; % make surface semi-transparent
s.EdgeColor = 'none'; % no grid lines
hold on
contour3(muA_mesh, c_mesh, hc_change, [0 0], 'k', 'LineWidth', 2)% Add black contour where hc_change = 0
offset = 0.12; % vertical offset above surface for better visibility
text(muA_mesh(4,1), c_mesh(4,1), hc_change(4,1) + offset, 'Employment Effect Dominates', ...
    'Color', 'black', 'FontSize', 7, 'FontWeight', 'bold', 'BackgroundColor', 'white', 'Margin', 2)
text(muA_mesh(4,3), c_mesh(4,3), hc_change(4,3) + offset, 'Investment Disincentive Effect Dominates', ...
    'Color', 'black', 'FontSize', 7, 'FontWeight', 'bold', 'BackgroundColor', 'white', 'Margin', 2)
xlabel('\mu_A Value')
xticks([0 2 4 6])
xticklabels({'0.1','0.2','0.3','0.4'})
yticks([0 2 4 6])
yticklabels({'0.6','0.5','0.4','0.3'})
ylabel('Vacancy Cost (c)')
zlabel('% Change Avg. Human Capital')
hold off
saveas(gcf,'Figure8','epsc')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Appendices %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
param = param_orig;  % Reset parameter values
%=======================================================================%
%----------------------- Appendix E1 -----------------------------------%
%=======================================================================%

indicator_policy=0; % Useful later when considering effects of policy changes
% Solve policy and value functions and steady-state distribution
[A_staffing,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,~,opt_p_N,opt_p_O,~,opt_x_N,opt_x_O,opt_eta_U,~,~,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          

sim_decomp_size = 50000;     % Number of individuals to follow only for one year in this case
[wp_decomp,wp_nc_decomp,wp_coll_decomp,wg_decomp,wg_nc_decomp,wg_coll_decomp] = DO_sim_decomp_func(dist_N,sim_decomp_size,param,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,opt_p_N,opt_p_O,opt_x_N,opt_x_O,opt_eta_U);

% Wage Penalty Results
wp_dec(1,1) = wp_decomp(1);
wp_dec(2,1) = wp_nc_decomp(1);
wp_dec(3,1) = wp_coll_decomp(1);

% Wage Growth Results
wg_dec(1,1) = wg_decomp(1);
wg_dec(2,1) = wg_nc_decomp(1);
wg_dec(3,1) = wg_coll_decomp(1);

%- Table 11: Wage Penalty - Comparison with Model Without Selection
fprintf('Wage Penalty Without Selection Effects\n');
T = table(wp_model,wp_dec,'RowName',{'Total','Only Non-College','Only College'},'VariableNames',{'Full Model','Model Without Selection'});
disp(T)

row_labels = {'Total', 'Non-College Only', 'College Only'};
fid = fopen('Table11.tex', 'w');
fprintf(fid, '\\begin{table}[h]\n');
fprintf(fid, '\\centering\n');
fprintf(fid, '\\caption{Wage Penalty: Full Model Versus Model Without Selection}\n');
fprintf(fid, '\\label{tab:wp_no_selection}\n');
fprintf(fid, '\\begin{tabular}{ c c c }\n');
fprintf(fid, '\\hline\n');
fprintf(fid, '    & Full Model (See Table 3) & Model Without Selection \\\\ \n');
fprintf(fid, '\\hline\n');
for i = 1:3
    fprintf(fid, ' %s & %.3f & %.3f \\\\  \n', row_labels{i}, wp_model(i), wp_dec(i));
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%- Table 12: Wage Growth - Comparison with Model Without Selection
fprintf('Wage Growth Without Selection Effects\n');
T = table(wg_model,wg_dec,'RowName',{'Total','Only Non-College','Only College'},'VariableNames',{'Full Model','Model Without Selection'});
disp(T)

row_labels = {'Total', 'Non-College Only', 'College Only'};
fid = fopen('Table12.tex', 'w');
fprintf(fid, '\\begin{table}[h]\n');
fprintf(fid, '\\centering\n');
fprintf(fid, '\\caption{Wage Growth Difference: Full Model Versus Model Without Selection}\n');
fprintf(fid, '\\label{tab:wg_no_selection}\n');
fprintf(fid, '\\begin{tabular}{ c c c }\n');
fprintf(fid, '\\hline\n');
fprintf(fid, '    & Full Model (See Table 4) & Model Without Selection \\\\ \n');
fprintf(fid, '\\hline\n');
for i = 1:3
    fprintf(fid, ' %s & %.3f & %.3f \\\\  \n', row_labels{i}, wg_model(i), wg_dec(i));
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);
%=======================================================================%
%----------------------- Appendix D5 -----------------------------------%
%=======================================================================%

% First Compute Steady State Comparison Values
npolicy = 2; % Number of policies considered
urate_ss = zeros(1,npolicy);
GDP_ss = zeros(1,npolicy);
avg_hc_ss = zeros(1,npolicy);
util_ss = zeros(1,npolicy);
pc_DO_ss = zeros(1,npolicy);

for policy = 1:npolicy
    
    if policy==2 % Need to re-solve the model in this case
        muA=1;
        param = [nu beta alpha d ell lambda chi b muA c delta h_sigma pr_phi z];
        % Solve policy and value functions and steady-state distribution
        [A_staffing,~,~,~,~,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,~,~,~,~,~,~,~,~,~,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          
    end
    
    % Compute and save moments of interest for comparison
    %-- Unemployment rate
    urate_ss(policy) = (sum(sum(dist_U(:,:)))/( sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:))) + sum(sum(dist_U(:,:))) ))*100;

    %----------- Output from Employment (GDP)
    for ind_phi = 1:nphi
        for ind_h = 1:nh
            h = h_grid(ind_h);
            GDP_ss(policy) = GDP_ss(policy) + (z(ind_phi)*h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
        end
    end

    %---------- Average Human Capital
    h_sum_total = 0;
    for ind_phi = 1:nphi
        for ind_h = 1:nh
            h = h_grid(ind_h);
            avg_hc_ss(policy) = avg_hc_ss(policy) + h*(dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
            h_sum_total = h_sum_total + dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h);
        end
    end
    avg_hc_ss(policy) = avg_hc_ss(policy)/h_sum_total;
   
    %--------- Total Utility
    util = 0;
    for ind_phi = 1:nphi
        for ind_h = 1:nh
            h = h_grid(ind_h);
            hp_N = h_grid(opt_ind_hp_N(ind_phi,ind_h));
            hp_O = h_grid(opt_ind_hp_O(ind_phi,ind_h));
            % Compute production minus investment cost
            util = util + dist_U(ind_phi,ind_h)*(b*z(ind_phi)*(h^alpha))...
                + z(ind_phi)*(h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h))...
                - (((hp_N - h*(1-d))*chi)^gamma)*dist_N(ind_phi,ind_h) - (((hp_O - h*(1-d))*chi)^gamma)*dist_O(ind_phi,ind_h);
            % Now subtract out vacancy posting costs
            if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3) % Search in unemployment for non-outsourced job
                if (~isnan(opt_theta_U(ind_phi,ind_h)))
                    util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*c;
                end
            else                                              % Search in unemployment for outsourced job
                if (~isnan(opt_theta_U(ind_phi,ind_h)))
                    util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*(c + A_staffing(ind_phi,ind_h));
                end
            end

            if (GN(ind_phi,ind_h)==1 || GN(ind_phi,ind_h)==3) % Search from non-outsourced job for non-outsourced job
                if (~isnan(opt_theta_N(ind_phi,ind_h)))
                    util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*c;
                end
            else                                              % Search from outsourced job for non-outsourced job
                if (~isnan(opt_theta_N(ind_phi,ind_h)))
                    util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                end
            end

            if (GO(ind_phi,ind_h)==1 || GO(ind_phi,ind_h)==3) % Search from outsourced job for non-outsourced job
                if (~isnan(opt_theta_O(ind_phi,ind_h)))
                    util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*c;
                end
            else
                if (~isnan(opt_theta_O(ind_phi,ind_h)))
                    util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                end
            end
        end
    end
    util_ss(policy) = util;

    % Percent Outsourced
    pc_DO_ss(policy) = (( sum(sum(dist_O(:,:))) )/(sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:)))))*100.0;
  
end % end loop over policy

% Now Compute Values along Transition Path
indicator_policy=0; % Useful later when considering effects of policy changes
% Solve policy and value functions and steady-state distribution
param = param_orig;
[A_staffing,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,opt_x_U,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy);          

param_tr = param_orig;
param_tr(9) = 1;       % Set muA=1 so there is no outsourcing
indicator_policy=0;
[util_tr,pc_DO_tr,urate_tr,avg_hc_tr,GDP_tr] = DO_transitions_func(dist_U,dist_N,dist_O,param_tr,indicator_policy);

ntr = 100;
t_grid=linspace(1,ntr,ntr);

pc_util_tr = zeros(1,ntr);
pc_util_ss = zeros(1,ntr);

pc_avgh_tr = zeros(1,ntr);
pc_avgh_ss = zeros(1,ntr);

pc_GDP_tr = zeros(1,ntr);
pc_GDP_ss = zeros(1,ntr);
for t = 1:ntr
    pc_util_tr(t) = ((util_tr(t) - util_ss(1))/util_ss(1))*100;
    pc_util_ss(t) = ((util_ss(2) - util_ss(1))/util_ss(1))*100;

    pc_avgh_tr(t) = ((avg_hc_tr(t) - avg_hc_ss(1))/avg_hc_ss(1))*100;
    pc_avgh_ss(t) = ((avg_hc_ss(2) - avg_hc_ss(1))/avg_hc_ss(1))*100;

    pc_GDP_tr(t) = ((GDP_tr(t) - GDP_ss(1))/GDP_ss(1))*100;
    pc_GDP_ss(t) = ((GDP_ss(2) - GDP_ss(1))/GDP_ss(1))*100;
end

figure
plot(t_grid,pc_util_tr,'b','LineWidth',2)
hold on
yline(0,'-.k','LineWidth',2)
plot(t_grid,pc_util_ss,'--r','LineWidth',2)
legend('Transition Path','Original Steady State','No Outsourcing Steady State')
xlabel('Quarters After Policy Change')
ylabel('Total Utility: % Change Relative to Original SS')
ylim([-0.2 0.01])
yticks([-0.2 -0.18 -0.16 -0.14 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0])
yticklabels({'-0.20%','-0.18%','-0.16%','-0.14%','-0.12%','-0.10%','-0.08%','-0.06%','-0.04%','-0.02%','0%'})
xticks([0 20 40 60 80 100])
saveas(gcf,'Figure9','epsc')

figure
plot(t_grid,pc_DO_tr,'b','LineWidth',2)
hold on
yline(pc_DO_ss(1),'-.k','LineWidth',2)
yline(pc_DO_ss(2),'--r','LineWidth',2)
legend('Transition Path','Original Steady State','No Outsourcing Steady State')
xlabel('Quarters After Policy Change')
ylabel('Percent of Employees in Outsourced Jobs')
yticks([0 0.5 1 1.5 2 2.5 3])
yticklabels({'0%','0.5%','1%','1.5%','2%','2.5%','3%'})
xticks([0 20 40 60 80 100])
saveas(gcf,'Figure10a','epsc')

figure
plot(t_grid,urate_tr,'b','LineWidth',2)
hold on
yline(urate_ss(1),'-.k','LineWidth',2)
yline(urate_ss(2),'--r','LineWidth',2)
legend('Transition Path','Original Steady State','No Outsourcing Steady State','Location','SouthEast')
xlabel('Quarters After Policy Change')
ylabel('Unemployment Rate')
ylim([4.85 5.2])
yticks([4.85 4.9 4.95 5 5.05 5.1 5.15 5.2])
yticklabels({'4.85%','4.9%','4.95%','5%','5.05%','5.1%','5.15%','5.2%'})
xticks([0 20 40 60 80 100])
saveas(gcf,'Figure10b','epsc')

figure
plot(t_grid,pc_avgh_tr,'b','LineWidth',2)
hold on
yline(0,'-.k','LineWidth',2)
plot(t_grid,pc_avgh_ss,'--r','LineWidth',2)
legend('Transition Path','Original Steady State','No Outsourcing Steady State','Location','NorthEast')
xlabel('Quarters After Policy Change')
ylabel('Avg. Human Capital: % Change Relative to Original SS')
ylim([-0.18 0.01])
yticks([-0.18 -0.16 -0.14 -0.12 -0.1 -0.08 -0.06 -0.04 -0.02 0])
yticklabels({'-0.18%','-0.16%','-0.14%','-0.12%','-0.1%','-0.08%','-0.06%','-0.04%','-0.02%','0%'})
xticks([0 20 40 60 80 100])
saveas(gcf,'Figure10c','epsc')

figure
plot(t_grid,pc_GDP_tr,'b','LineWidth',2)
hold on
yline(0,'-.k','LineWidth',2)
plot(t_grid,pc_GDP_ss,'--r','LineWidth',2)
legend('Transition Path','Original Steady State','No Outsourcing Steady State','Location','NorthEast')
xlabel('Quarters After Policy Change')
ylabel('GDP: % Change Relative to Original SS')
ylim([-0.25 0.01])
yticks([-0.25 -0.2 -0.15 -0.1 -0.05 0])
yticklabels({'-0.25%','-0.2%','-0.15%','-0.1%','-0.05%','0%'})
xticks([0 20 40 60 80 100])
saveas(gcf,'Figure10d','epsc')

%=======================================================================%
%----------------------- Appendix E3 -----------------------------------%
%=======================================================================%
indicator_policy=0;
c_params = param_calibrated';
eps_machine = eps;
n_params = length(c_params);
J = zeros(n_params, n_params);

% Consider 1% increase or decrease in param values 
increment = zeros(1,n_params);
for i = 1:n_params
    increment(i) = param_calibrated(i)*0.01;
end

mm_inc = zeros(n_params,n_params);
mm_dec = zeros(n_params,n_params);
pc_mm_inc = zeros(n_params,n_params);
pc_mm_dec = zeros(n_params,n_params);
for i=1:n_params
    
    dp = zeros(n_params, 1);
    dp(i) = increment(i);

    %--- Compute moments with inputs c_params + dp ---%
    temp_c_params = c_params + dp;
    param_temp = [param(1) param(2) temp_c_params(8) temp_c_params(6) param(5) temp_c_params(1) temp_c_params(7) temp_c_params(2) temp_c_params(3) temp_c_params(4) temp_c_params(5) temp_c_params(10) param(13) param(14) param(15) temp_c_params(9)];

    % Solve policy and value functions and steady-state distribution
    [~,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,~,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O,~,~,~,dist_U,dist_N,dist_O] = DO_main_func(param_temp,indicator_policy);          

    % Simulate wages 
    simul_size = 20000;     % Number of individuals to follow in simulated sample
    [simul,~] = DO_sim_func(simul_size,param_temp,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O);

    % Generate Model Moments
    [~,~,~,~,~,~,~,~,~,~,mm,~,~,~,~,~,~] = DO_moments_func(h_grid,simul,simul_size,param_temp,dist_U,dist_N,dist_O,opt_ind_hp_N,opt_ind_hp_O,opt_p_U,opt_p_N,opt_p_O);

    mm_inc(:,i) = mm;

    %--- Compute moments with inputs c_params - dp
    temp_c_params = c_params - dp;
    param_temp = [param(1) param(2) temp_c_params(8) temp_c_params(6) param(5) temp_c_params(1) temp_c_params(7) temp_c_params(2) temp_c_params(3) temp_c_params(4) temp_c_params(5) temp_c_params(10) param(13) param(14) param(15) temp_c_params(9)];

    % Solve policy and value functions and steady-state distribution
    [A_staffing,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,opt_x_U,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param_temp,indicator_policy);          

    % Simulate wages 
    simul_size = 20000;     % Number of individuals to follow in simulated sample
    [simul,~] = DO_sim_func(simul_size,param_temp,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O);

    % Generate Model Moments
    [~,~,wage_sample_N1,wage_sample_N2,wage_sample_O1,wage_sample_O2,pr,pr_qrtr,pr_noncollege,pr_college,mm,wp_results,wp_noncollege_results,wp_college_results,wg_results,wg_noncollege,wg_college] = DO_moments_func(h_grid,simul,simul_size,param_temp,dist_U,dist_N,dist_O,opt_ind_hp_N,opt_ind_hp_O,opt_p_U,opt_p_N,opt_p_O);

    mm_dec(:,i) = mm;

    for j = 1:n_params
        pc_mm_inc(j,i) = ((mm_inc(j,i) - mm_orig(j))/abs(mm_orig(j)))*100;

        pc_mm_dec(j,i) = ((mm_dec(j,i) - mm_orig(j))/abs(mm_orig(j)))*100;
    end

end


rowNames = ["Job-to-job transition rate", "Unemployment rate (\%)", "Percent outsourced", ...
    "(Outsourced w/out)/(w/ Bachelor's)", "Average quarterly EU rate", ...
    "Avg. wage loss after unemp. spell (\%)", "Wage dispersion: p90/p50", "Avg. annual real wage growth (\%)", ...
    "College wage premium", "Age 20 Unemployment rate (\%)"];
colNames = ["$\lambda$", "$b_c$", "$\mu^A$", "$c$", "$\delta$","$d$", ...
            "$\chi$", "$\alpha$", "$z(\phi_H)$", "$\sigma_h$"];

fileID = fopen('Table15.tex', 'w');
fprintf(fileID, '\\begin{table}[h]\n\\centering\n');
fprintf(fileID, '\\caption{Percent Change in Moments from 1 Percent Parameter Increases}\n');
fprintf(fileID, '\\label{tab:pc_mm_table}\n');
fprintf(fileID, '\\scalebox{0.72}{\n\\begin{tabular}{l%s}\n', repmat('c', 1, 10));
fprintf(fileID, '\\hline\n');
fprintf(fileID, ' & ');
fprintf(fileID, '%s', colNames(1));
for j = 2:length(colNames)
    fprintf(fileID, ' & %s', colNames(j));
end
fprintf(fileID, ' \\\\\n\\hline\n');
for i = 1:10
    fprintf(fileID, '%s', rowNames(i));
    for j = 1:10
        fprintf(fileID, ' & %.3f', pc_mm_inc(i, j)); 
    end
    fprintf(fileID, ' \\\\\n');
end

fprintf(fileID, '\\hline\n\\end{tabular}}\n\\end{table}\n');
fclose(fileID);

fileID = fopen('Table16.tex', 'w');
fprintf(fileID, '\\begin{table}[h]\n\\centering\n');
fprintf(fileID, '\\caption{Percent Change in Moments from 1 Percent Parameter Decreases}\n');
fprintf(fileID, '\\label{tab:pc_mm_dec_table}\n');
fprintf(fileID, '\\scalebox{0.72}{\n\\begin{tabular}{l%s}\n', repmat('c', 1, 10));
fprintf(fileID, '\\hline\n');
fprintf(fileID, ' & ');
fprintf(fileID, '%s', colNames(1));
for j = 2:length(colNames)
    fprintf(fileID, ' & %s', colNames(j));
end
fprintf(fileID, ' \\\\\n\\hline\n');
for i = 1:10
    fprintf(fileID, '%s', rowNames(i));
    for j = 1:10
        fprintf(fileID, ' & %.3f', pc_mm_dec(i, j));  
    end
    fprintf(fileID, ' \\\\\n');
end

fprintf(fileID, '\\hline\n\\end{tabular}}\n\\end{table}\n');
fclose(fileID);

%=======================================================================%
%----------------------- Appendix E1 -----------------------------------%
%=======================================================================%
param = param_orig;
[mm_E1] = DO_Appendix_E1(param);
%--------%
rowNames = {'$\lambda$','$b_c$','$\zeta^F$','$c$','$\delta$','$d$','$\chi$','$\alpha$', ...
            '$z(\phi_H)$','$\sigma_h$'};
fprintf('Jointly Calibrated Parameters\n');
targets = ["Job-to-job transition rate";"Unemployment rate (\%)";"Percent outsourced via staffing agency";"Ratio outsourced w/out to outsourced w/ Bachelor's";"Average quarterly EU rate";"Avg. wage loss after unemployment spell (\%)";...
    "Wage dispersion: p90/p50";"Avg. annual real wage growth (\%)";"College wage premium";"Age 20 Unemployment rate (\%)"];
T = table(param_calibrated',targets,mm_baseline,mm_E1,'VariableNames',{'Estimate','Targeted Moment','Main_Text','E1'},'RowName',rowNames); 
disp(T)     

% ---- Save Table 17
fid = fopen('Table17.tex', 'w');
fprintf(fid, '\\begin{table}[htbp]\n\\centering\n');
fprintf(fid, '\\caption{Jointly Calibrated Parameters}\n');
fprintf(fid, '\\scalebox{0.82}{%%\n');
fprintf(fid, '\\begin{tabular}{l c c l c c}\n');
fprintf(fid, '\\hline\n');
fprintf(fid, ' & Main Text Estimate & E.1 Estimate & Targeted Moment & Main Text & E.1 \\\\\n');
fprintf(fid, '\\hline\n');
for i = 1:height(T)
    rowName    = T.Properties.RowNames{i};    % Already contains LaTeX ($...$)
    est_main   = T.Estimate(i);     % numeric
    est_e1     = T.Estimate(i);            % numeric
    moment     = T.("Targeted Moment"){i};    % string
    data_main  = sprintf('%.3f', T.Main_Text(i));
    data_e1    = sprintf('%.3f', T.E1(i));

    fprintf(fid, '%s & %.6f & %.6f & %s & %s & %s \\\\\n', ...
        rowName, est_main, est_e1, moment, data_main, data_e1);
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}}\n');
fprintf(fid, '\\label{tab:D1_calibration}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%=======================================================================%
%----------------------- Appendix E2 -----------------------------------%
%=======================================================================%

% Table 18 - Table 19
%---- Model is re-calibrated ----%
pr_phi(1) = 1.0-0.3084251;    % Probability of entering the model with phi=1 (low fixed productivity type)
pr_phi(2) = 0.3084251;        % Probability of entering the model with phi=2 (high fixed productivity type)

z(1) = 1;                     % Productivity parameter f(phi,h) = z(phi)h^(alpha)
z(2) = 1.171936;

param_E2 = [0.00622206 0.99*(1-nu) 0.61031 0.028695 1.27 0.423303 4.043024 0.461709 0.110444 0.684773 0.114167 0.239816 pr_phi z];
param_E2_est = [0.423303 0.461709 0.110444 0.684773 0.114167 0.028695 4.043024 0.61031 1.171936 0.239816];
psi = 0.05;
[mm,opt_hp_N,cdf_N,cdf_O,avg_diff_N,avg_diff_O,policy_change_table] = DO_Appendix_E2(param_E2,psi);


%---- Table 18
rowNames = {'$\lambda$','$b_c$','$\mu^A$','$c$','$\delta$','$d$','$\chi$','$\alpha$', ...
            '$z(\phi_H)$','$\sigma_h$'};
fprintf('Jointly Calibrated Parameters\n');
targets = ["Job-to-job transition rate";"Unemployment rate (\%)";"Percent outsourced via staffing agency";"Ratio outsourced w/out to outsourced w/ Bachelor's";"Average quarterly EU rate";"Avg. wage loss after unemployment spell (\%)";...
    "Wage dispersion: p90/p50";"Avg. annual real wage growth (\%)";"College wage premium";"Age 20 Unemployment rate (\%)"];
T = table(param_E2_est',targets,dm,mm,'VariableNames',{'Estimate','Targeted Moment','Data','Model'},'RowName',rowNames); 
disp(T)     

% Open file for writing LaTeX table
fid = fopen('Table18.tex', 'w');
fprintf(fid, '\\begin{table}[htbp]\n\\centering\n');
fprintf(fid, '\\caption{Jointly Calibrated Parameters - With Job-Specific Human Capital ($\\psi=0.5$)}\n');
fprintf(fid, '\\scalebox{0.92}{%%\n');
fprintf(fid, '\\begin{tabular}{l c l c c}\n');
fprintf(fid, '\\hline\n');
fprintf(fid, 'Parameter & Estimate & Targeted Moment & Data & Model \\\\\n');
fprintf(fid, '\\hline\n');
% Loop through table rows
for i = 1:height(T)
    rowName  = T.Properties.RowNames{i};                 % Already contains LaTeX ($...$)
    estimate = T.Estimate(i);                            % numeric
    moment   = T.("Targeted Moment"){i};                 % string
    dataStr  = sprintf('%.3f', T.Data(i));               % format Data
    modelStr = sprintf('%.3f', T.Model(i));              % format Model
    fprintf(fid, '%s & %.6f & %s & %s & %s \\\\\n', ...
        rowName, estimate, moment, dataStr, modelStr);
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}}\n');
fprintf(fid, '\\label{Appendix_psi_5pc}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%---- Figure 11a            
figure
plot(h_grid,opt_hp_N(1,:),'LineWidth',2)
hold on
plot(h_grid,opt_hp_N(2,:),'Color',[0.9290 0.6940 0.1250],'LineWidth',2)
plot(h_grid,h_grid,':k','LineWidth',1)
xlim([0.0 7])
yticks([0 1 2 3 4 5 6 7])
legend('Non-College','College','45-Degree Line','Location','NorthWest')
xlabel('Current Human Capital h (Relative to New Entrant Mean Normalized to 1)')
ylabel('Optimal Next Period Human Capital')
saveas(gcf,'Figure11a','epsc')   

%---- Figure 11b  
figure
plot(h_grid,cdf_N(1,:),'Color',[0 0.4470 0.7410],'LineWidth',2)
hold on
plot(h_grid,cdf_O(1,:),'Color',[0.8500 0.3250 0.0980],'LineWidth',2)
plot(h_grid,cdf_N(2,:),'Color',[0.9290 0.6940 0.1250],'LineWidth',2)
plot(h_grid,cdf_O(2,:),'Color',[0.4940 0.1840 0.5560],'LineWidth',2)
legend('Non-College: Not Outsourced','Non-College: Outsourced','College: Not Outsourced','College: Outsourced','location','NorthWest')
ylabel('CDF')
xlabel('Human Capital (Relative to New Entrant Mean Normalized to 1)')
ylim([0 1])
xlim([0.0 7])
saveas(gcf,'Figure11b','epsc')  

%---- Figure 12
h_groups = linspace(1,length(avg_diff_N(1,:)),length(avg_diff_N(1,:)));
figure
scatter(h_groups,avg_diff_N(1,:),'b','filled')
hold on
plot(h_groups,avg_diff_N(1,:),':b','HandleVisibility','off')
scatter(h_groups,avg_diff_O(1,:),'s','r','filled')
plot(h_groups,avg_diff_O(1,:),':r','HandleVisibility','off')
scatter(h_groups,avg_diff_N(2,:),'k','filled')
plot(h_groups,avg_diff_N(2,:),':k','HandleVisibility','off')
scatter(h_groups,avg_diff_O(2,:),'s','MarkerEdgeColor',[0.9290 0.6940 0.1250],'MarkerFaceColor',[0.9290 0.6940 0.1250])
plot(h_groups,avg_diff_O(2,:),'Color',[0.9290 0.6940 0.1250],'LineStyle',':','HandleVisibility','off')
yline(0,'k')
legend('Non-College: Not Outsourced','Non-College: Outsourced','College: Not Outsourced','College: Outsourced')
ylabel('Optimal Human Capital Increase: h`-h')
xlabel('Current Human Capital (h)')
xticks([1 3 5 7 9 11 ])
xticklabels({'[0,0.25)','[0.5,0.75)','[1,1.25)','[1.5,1.75)','[2,2.25)','[2.5,2.75)'})
ylim([0.1 0.305])
xlim([1 12])
box on
saveas(gcf,'Figure12','epsc') 

%---- Table 19
fprintf('Table 19: Steady State Comparison: Aggregate Effects of Removing Staffing Agencies\n');
T = table(total_results,policy_change_table,'RowName',{'Unemployment','Output','Average_HC','Average_Wage','Total_Utility','Entrant_Value'},'VariableNames',{'Main Text','With Job-Specific Human Capital'});
disp(T)

% Write table 19 to a LaTeX file
row_labels = {'Unemployment rate';'Output from employment';'Average human capital';
    'Average wage';'Total consumption';'Entrant value'};

% Units for each row (e.g., "pp" or "\%%" for LaTeX)
row_units = {'pp';'\%';'\%';'\%';'\%';'\%'};

fileID = fopen('Table19.tex','w');
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Steady State Comparison: Aggregate Effects of Removing Staffing Agencies ($\\psi=0.05$)}\n');
fprintf(fileID, '\\label{tab:Appendix_5pc_counterfactual}\n');
fprintf(fileID, '\\begin{tabular}{ c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '     & Main Text (Table 7) & With Job-Specific Human Capital \\\\ \n');
fprintf(fileID, '\\hline\n');

for i = 1:6
    % Print with sign (+/-), 3 decimal places, and units
    fprintf(fileID, ' %s & %+.3f%s & %+.3f%s \\\\\n', ...
        row_labels{i}, ...
        total_results(i), row_units{i}, ...
        policy_change_table(i), row_units{i});
end
fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);

% Table 20 - Table 21
%---- Model is re-calibrated ----%
z(1) = 1;                     % Productivity parameter f(phi,h) = z(phi)h^(alpha)
z(2) = 1.144724;
param_E2 = [0.00622206 0.99*(1-nu) 0.612048 0.006835 1.27 0.400353 4.065891 0.519997 0.111744 0.712333 0.120635 0.23514 pr_phi z];
param_E2_est = [0.400353 0.519997 0.111744 0.712333 0.120635 0.006835 4.065891 0.612048 1.144724 0.23514];
psi = 0.1;
[mm,opt_hp_N,cdf_N,cdf_O,~,~,policy_change_table] = DO_Appendix_E2(param_E2,psi);

%---- Table 20
rowNames = {'$\lambda$','$b_c$','$\mu^A$','$c$','$\delta$','$d$','$\chi$','$\alpha$', ...
            '$z(\phi_H)$','$\sigma_h$'};
fprintf('Jointly Calibrated Parameters\n');
targets = ["Job-to-job transition rate";"Unemployment rate (\%)";"Percent outsourced via staffing agency";"Ratio outsourced w/out to outsourced w/ Bachelor's";"Average quarterly EU rate";"Avg. wage loss after unemployment spell (\%)";...
    "Wage dispersion: p90/p50";"Avg. annual real wage growth (\%)";"College wage premium";"Age 20 Unemployment rate (\%)"];
T = table(param_E2_est',targets,dm,mm,'VariableNames',{'Estimate','Targeted Moment','Data','Model'},'RowName',rowNames); 
disp(T)     

% Open file for writing LaTeX table
fid = fopen('Table20.tex', 'w');
fprintf(fid, '\\begin{table}[htbp]\n\\centering\n');
fprintf(fid, '\\caption{Jointly Calibrated Parameters - With Job-Specific Human Capital ($\\psi=0.1$)}\n');
fprintf(fid, '\\scalebox{0.92}{%%\n');
fprintf(fid, '\\begin{tabular}{l c l c c}\n');
fprintf(fid, '\\hline\n');
fprintf(fid, 'Parameter & Estimate & Targeted Moment & Data & Model \\\\\n');
fprintf(fid, '\\hline\n');

% Loop through table rows
for i = 1:height(T)
    rowName  = T.Properties.RowNames{i};                 % Already contains LaTeX ($...$)
    estimate = T.Estimate(i);                            % numeric
    moment   = T.("Targeted Moment"){i};                 % string
    dataStr  = sprintf('%.3f', T.Data(i));               % format Data
    modelStr = sprintf('%.3f', T.Model(i));              % format Model
    fprintf(fid, '%s & %.6f & %s & %s & %s \\\\\n', ...
        rowName, estimate, moment, dataStr, modelStr);
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}}\n');
fprintf(fid, '\\label{Appendix_D_pc10:calibration}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%---- Table 21
fprintf('Table 21: Steady State Comparison: Aggregate Effects of Removing Staffing Agencies\n');
T = table(total_results,policy_change_table,'RowName',{'Unemployment','Output','Average_HC','Average_Wage','Total_Utility','Entrant_Value'},'VariableNames',{'Main Text','With Job-Specific Human Capital'});
disp(T)

% Write table 21 to a LaTeX file
row_labels = {'Unemployment rate';'Output from employment';'Average human capital';
    'Average wage';'Total consumption';'Entrant value'};

% Units for each row (e.g., "pp" or "\%%" for LaTeX)
row_units = {'pp';'\%';'\%';'\%';'\%';'\%'};

fileID = fopen('Table21.tex','w');
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Steady State Comparison: Effects of Removing Staffing Agencies ($\\psi=0.1$)}\n');
fprintf(fileID, '\\label{tab:Appendix_10pc_hc}\n');
fprintf(fileID, '\\begin{tabular}{ c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '     & Main Text (Table 7)  & With Job-Specific Human Capital \\\\ \n');
fprintf(fileID, '\\hline\n');

for i = 1:length(row_labels)
    % Print with sign (+/-), 3 decimal places, and units
    fprintf(fileID, ' %s & %+.3f%s & %+.3f%s \\\\\n', ...
        row_labels{i}, ...
        total_results(i), row_units{i}, ...
        policy_change_table(i), row_units{i});
end
fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);

%=======================================================================%
%----------------------- Appendix E3 -----------------------------------%
%=======================================================================%

[mm_E3,param_E3,policy_E3,perc_college_O,avg_h_O,separation_O,perc_college_N,avg_h_N,separation_N,avg_diff_N,avg_diff_O,avg_diff_L,avg_diff_H] = DO_Appendix_E3;

%---- Table 22 ----%
dm_E3 = [0.066;4.835;2.618;1.366;4.650;0.034;3.122;-4.400;2.133;2.357;1.791;11.779];
rowNames = {'$\lambda$','$b_c$','$\mu^A$','$c_L$','$c_H$','$\delta_L$','$\delta_H$','$d$','$\chi$','$\alpha$', ...
    '$z(\phi_H)$','$\sigma_h$'};
fprintf('Table 20: Jointly Calibrated Parameters\n');
targets = ["Job-to-job transition rate";"Unemployment rate (\%)";"Percent outsourced via staffing agency";"Ratio outsourced w/out to outsourced w/ Bachelor's";...
    "Job vacancy rate";"Average quarterly EU rate";"Avg. Outsourced/non-outsourced separation rate";"Avg. wage loss after unemployment spell (\%)";...
    "Wage dispersion: p90/p50";"Avg. annual real wage growth (\%)";"College wage premium";"Age 20 Unemployment rate (\%)"];
T = table(param_E3,targets,dm_E3,mm_E3,'VariableNames',{'Estimate','Targeted Moment','Data','Model'},'RowName',rowNames); 
disp(T) 

fid = fopen('Table22.tex', 'w');
fprintf(fid, '\\begin{table}[h!]\n\\centering\n');
fprintf(fid, '\\caption{Jointly Calibrated Parameters - Multiple Job Separation Types}\n');
fprintf(fid, '\\scalebox{0.85}{\n');
fprintf(fid, '\\begin{tabular}{l c l c c}\n');
fprintf(fid, '\\hline\n');
fprintf(fid, 'Parameter & Estimate & Targeted Moment & Data & Model \\\\\n');
fprintf(fid, '\\hline\n');

% Loop through table rows
for i = 1:height(T)
    rowName  = T.Properties.RowNames{i};                 % Already contains LaTeX ($...$)
    estimate = T.Estimate(i);                            % numeric
    moment = char(T.("Targeted Moment")(i));
    dataStr  = sprintf('%.3f', T.Data(i));               % format Data
    modelStr = sprintf('%.3f', T.Model(i));              % format Model
    fprintf(fid, '%s & %.6f & %s & %s & %s \\\\\n', ...
        rowName, estimate, moment, dataStr, modelStr);
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}}\n');
fprintf(fid, '\\label{Appendix_D3:calibration}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%---- Table 23 ----%
fprintf('Table 23: Selection Into Outsourcing: Model with Mulitple Job Separation Types\n');
T = table([perc_college_O;avg_h_O;separation_O],[perc_college_N;avg_h_N;separation_N],'VariableNames',{'Outsourced','Non_Outsourced'},'RowName',{'Percent_with_Bachelors','Average_h','Separation_Rate'});
disp(T)
     
% Open file to write LaTeX table
fid = fopen('Table23.tex', 'w');
% Start LaTeX table
fprintf(fid, '\\begin{table}[h]\n');
fprintf(fid, '\\centering\n');
fprintf(fid, '\\caption{Selection Into Outsourcing: Model with Multiple Job Separation Types}\n');
fprintf(fid, '\\label{tab:D3_Selection}\n');
fprintf(fid, '\\scalebox{0.9}{\n');
fprintf(fid, '\\begin{tabular}{ c c c }\n');
fprintf(fid, '\\hline\n');
fprintf(fid, ' & Outsourced & Non-Outsourced \\\\\n');
fprintf(fid, '\\hline\n');

% Get row labels and format rows
rowLabels = T.Properties.RowNames;
for i = 1:height(T)
    rowName = rowLabels{i};
    
    % Custom LaTeX label replacements
    switch rowName
        case 'Percent_with_Bachelors'
            label = 'Percent with Bachelor''s ($\phi=H$)';
            outsourced = sprintf('%.3f\\%%', T.Outsourced(i));  % percent format
            non_outsourced = sprintf('%.3f\\%%', T.Non_Outsourced(i));
        case 'Average_h'
            label = 'Average h (Relative to New Entrant Average = 1)';
            outsourced = sprintf('%.3f', T.Outsourced(i));
            non_outsourced = sprintf('%.3f', T.Non_Outsourced(i));
        case 'Separation_Rate'
            label = 'Average quarterly EU rate';
            outsourced = sprintf('%.3f', T.Outsourced(i));
            non_outsourced = sprintf('%.3f', T.Non_Outsourced(i));
        otherwise
            label = rowName; % fallback
            outsourced = sprintf('%.3f', T.Outsourced(i));
            non_outsourced = sprintf('%.3f', T.Non_Outsourced(i));
    end
    
    % Write table row
    fprintf(fid, ' %s & %s & %s \\\\\n', label, outsourced, non_outsourced);
end

% Close LaTeX table
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%---- Figure 13 ----%
h_groups = linspace(1,length(avg_diff_N),length(avg_diff_N));
        
figure
scatter(h_groups,avg_diff_N,'b','filled')
hold on
plot(h_groups,avg_diff_N,':b','HandleVisibility','off')
scatter(h_groups,avg_diff_O,'s','r','filled')
plot(h_groups,avg_diff_O,':r','HandleVisibility','off')
yline(0,'k')
legend('Not Outsourced','Outsourced')
ylabel('Optimal Human Capital Increase: h`-h')
xlabel('Current Human Capital (h)')
xticks([1 3 5 7 9 11 ])
xticklabels({'[0,0.25)','[0.5,0.75)','[1,1.25)','[1.5,1.75)','[2,2.25)','[2.5,2.75)'})
axis tight
ylim([0.05 0.25])
box on
saveas(gcf,'Figure13a','epsc') 
        
figure
scatter(h_groups,avg_diff_L,'b','filled')
hold on
plot(h_groups,avg_diff_L,':b','HandleVisibility','off')
scatter(h_groups,avg_diff_H,'s','g','filled')
plot(h_groups,avg_diff_H,':g','HandleVisibility','off')
yline(0,'k')
legend('Low Separation Rate Job','High Separation Rate Job')
ylabel('Optimal Human Capital Increase: h`-h')
xlabel('Current Human Capital (h)')
xticks([1 3 5 7 9 11 ])
xticklabels({'[0,0.25)','[0.5,0.75)','[1,1.25)','[1.5,1.75)','[2,2.25)','[2.5,2.75)'})
axis tight
ylim([0.05 0.25])
box on
saveas(gcf,'Figure13b','epsc') 

%---- Table 24
fprintf('Table 24: Steady State Comparison: Effects of Removing Staffing Agencies\n');
T = table(total_results,policy_E3,'RowName',{'Unemployment','Output','Average_HC','Average_Wage','Total_Utility','Entrant_Value'},'VariableNames',{'Main Text','With Job-Specific Human Capital'});
disp(T)

% Write table 24 to a LaTeX file
row_labels = {'Unemployment rate';'Output from employment';'Average human capital';
    'Average wage';'Total consumption';'Entrant value'};

% Units for each row (e.g., "pp" or "\%%" for LaTeX)
row_units = {'pp';'\%';'\%';'\%';'\%';'\%'};

fileID = fopen('Table24.tex','w');
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Steady State Comparison: Aggregate Effects of Removing Staffing Agencies}\n');
fprintf(fileID, '\\label{Appendix_deltas_counterfactual}\n');
fprintf(fileID, '\\begin{tabular}{ c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '     & Main Text (Table 7) & Multiple Job Separation Types \\\\ \n');
fprintf(fileID, '\\hline\n');

for i = 1:6
    % Print with sign (+/-), 3 decimal places, and units
    fprintf(fileID, ' %s & %+.3f%s & %+.3f%s \\\\\n', ...
        row_labels{i}, ...
        total_results(i), row_units{i}, ...
        policy_E3(i), row_units{i});
end
fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);

%=======================================================================%
%----------------------- Appendix E4 -----------------------------------%
%=======================================================================%

[mm_E4,param_E4,policy_E4,perc_college_O,avg_h_O,alpha_O,perc_college_N,avg_h_N,alpha_N,avg_diff_N,avg_diff_O,avg_diff_L,avg_diff_H] = DO_Appendix_E4;
    
%---- Table 25 ----%
dm_E4 = [0.066;4.835;2.618;1.366;4.650;0.034;-4.400;2.133;2.357;0.874;1.791;11.779];
rowNames = {'$\lambda$','$b_c$','$\mu^A$','$c_L$','$c_H$','$\delta$','$d$','$\chi$','$\alpha_L$','$\alpha_H$', ...
    '$z(\phi_H)$','$\sigma_h$'};
fprintf('Table 23: Jointly Calibrated Parameters\n');
targets = ["Job-to-job transition rate";"Unemployment rate (\%)";"Percent outsourced via staffing agency";"Ratio outsourced w/out to outsourced w/ Bachelor's";...
    "Job vacancy rate";"Average quarterly EU rate";"Avg. wage loss after unemployment spell (\%)";...
    "Wage dispersion: p90/p50";"Avg. annual real wage growth (\%)";"Wage ratio: 25-34/35-25 to $<$25/25-34";...
    "College wage premium";"Age 20 Unemployment rate (\%)"];
T = table(param_E4,targets,dm_E4,mm_E4,'VariableNames',{'Estimate','Targeted Moment','Data','Model'},'RowName',rowNames); 
disp(T) 


fid = fopen('Table25.tex', 'w');
fprintf(fid, '\\begin{table}[h]\n\\centering\n');
fprintf(fid, '\\caption{Jointly Calibrated Parameters - Multiple Sectors}\n');
fprintf(fid, '\\scalebox{0.9}{\n');
fprintf(fid, '\\begin{tabular}{l c l c c}\n');
fprintf(fid, '\\hline\n');
fprintf(fid, 'Parameter & Estimate & Targeted Moment & Data & Model \\\\\n');
fprintf(fid, '\\hline\n');

% Loop through table rows
for i = 1:height(T)
    rowName  = T.Properties.RowNames{i};                 % Already contains LaTeX ($...$)
    estimate = T.Estimate(i);                            % numeric
    moment = char(T.("Targeted Moment")(i));
    dataStr  = sprintf('%.3f', T.Data(i));               % format Data
    modelStr = sprintf('%.3f', T.Model(i));              % format Model
    fprintf(fid, '%s & %.6f & %s & %s & %s \\\\\n', ...
        rowName, estimate, moment, dataStr, modelStr);
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}}\n');
fprintf(fid, '\\label{Appendix_D4:calibration}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);


%---- Table 26 ----%                    
fprintf('Table 26: Model Selection Into Outsourcing\n');
T = table([perc_college_O;avg_h_O;alpha_O],[perc_college_N;avg_h_N;alpha_N],'VariableNames',{'Outsourced','Non_Outsourced'},'RowName',{'Percent_with_Bachelors','Average_h','Sectoral_alpha'});
disp(T)
     
% Open file to write LaTeX table
fid = fopen('Table26.tex', 'w');
% Start LaTeX table
fprintf(fid, '\\begin{table}[h]\n');
fprintf(fid, '\\centering\n');
fprintf(fid, '\\caption{Selection Into Outsourcing: Multiple Sectors}\n');
fprintf(fid, '\\label{tab:D4_Selection}\n');
fprintf(fid, '\\scalebox{0.9}{\n');
fprintf(fid, '\\begin{tabular}{ c c c }\n');
fprintf(fid, '\\hline\n');
fprintf(fid, ' & Outsourced & Non-Outsourced \\\\\n');
fprintf(fid, '\\hline\n');

rowLabels = T.Properties.RowNames;
for i = 1:height(T)
    rowName = rowLabels{i};
    
    % Custom LaTeX label replacements
    switch rowName
        case 'Percent_with_Bachelors'
            label = 'Percent with Bachelor''s ($\phi=H$)';
            outsourced = sprintf('%.3f\\%%', T.Outsourced(i));  % percent format
            non_outsourced = sprintf('%.3f\\%%', T.Non_Outsourced(i));
        case 'Average_h'
            label = 'Average h (Relative to New Entrant Average = 1)';
            outsourced = sprintf('%.3f', T.Outsourced(i));
            non_outsourced = sprintf('%.3f', T.Non_Outsourced(i));
        case 'Sectoral_alpha'
            label = 'Average Sectoral $\alpha$';
            outsourced = sprintf('%.3f', T.Outsourced(i));
            non_outsourced = sprintf('%.3f', T.Non_Outsourced(i));
        otherwise
            label = rowName; % fallback
            outsourced = sprintf('%.3f', T.Outsourced(i));
            non_outsourced = sprintf('%.3f', T.Non_Outsourced(i));
    end
    % Write table row
    fprintf(fid, ' %s & %s & %s \\\\\n', label, outsourced, non_outsourced);
end
fprintf(fid, '\\hline\n');
fprintf(fid, '\\end{tabular}}\n');
fprintf(fid, '\\end{table}\n');
fclose(fid);

%---- Figure 14 ----% 
h_groups = linspace(1,length(avg_diff_N),length(avg_diff_N));
                
figure
scatter(h_groups,avg_diff_N,'b','filled')
hold on
plot(h_groups,avg_diff_N,':b','HandleVisibility','off')
scatter(h_groups,avg_diff_O,'s','r','filled')
plot(h_groups,avg_diff_O,':r','HandleVisibility','off')
yline(0,'k')
legend('Not Outsourced','Outsourced')
ylabel('Optimal Human Capital Increase: h`-h')
xlabel('Current Human Capital (h)')
xticks([1 3 5 7 9 11 ])
xticklabels({'[0,0.25)','[0.5,0.75)','[1,1.25)','[1.5,1.75)','[2,2.25)','[2.5,2.75)'})
axis tight
ylim([0.05 0.25])
box on
saveas(gcf,'Figure14a','epsc') 

figure
scatter(h_groups,avg_diff_L,'b','filled')
hold on
plot(h_groups,avg_diff_L,':b','HandleVisibility','off')
scatter(h_groups,avg_diff_H,'s','g','filled')
plot(h_groups,avg_diff_H,':g','HandleVisibility','off')
yline(0,'k')
legend('Low \alpha','High \alpha')
ylabel('Optimal Human Capital Increase: h`-h')
xlabel('Current Human Capital (h)')
xticks([1 3 5 7 9 11 ])
xticklabels({'[0,0.25)','[0.5,0.75)','[1,1.25)','[1.5,1.75)','[2,2.25)','[2.5,2.75)'})
axis tight
ylim([0.05 0.25])
box on
saveas(gcf,'Figure14b','epsc')

%---- Table 27 ----% 
fprintf('Table 27: Steady State Comparison: Effects of Removing Staffing Agencies\n');
T = table(total_results,policy_E4,'RowName',{'Unemployment','Output','Average_HC','Average_Wage','Total_Utility','Entrant_Value'},'VariableNames',{'Main Text','With Job-Specific Human Capital'});
disp(T)

% Write table 27 to a LaTeX file
row_labels = {'Unemployment rate';'Output from employment';'Average human capital';'Average wage';
    'Total consumption';'Entrant value'};

% Units for each row (e.g., "pp" or "\%%" for LaTeX)
row_units = {'pp';'\%';'\%';'\%';'\%';'\%'};

fileID = fopen('Table27.tex','w');
fprintf(fileID, '\\begin{table}[h]\n');
fprintf(fileID, '\\centering\n');
fprintf(fileID, '\\caption{Steady State Comparison: Aggregate Effects of Removing Staffing Agencies}\n');
fprintf(fileID, '\\label{tab:D4_counterfactual}\n');
fprintf(fileID, '\\begin{tabular}{ c c c }\n');
fprintf(fileID, '\\hline\n');
fprintf(fileID, '     & Main Text (Table 7) & Multiple Sectors \\\\ \n');
fprintf(fileID, '\\hline\n');

for i = 1:6
    % Print with sign (+/-), 3 decimal places, and units
    fprintf(fileID, ' %s & %+.3f%s & %+.3f%s \\\\\n', ...
        row_labels{i}, ...
        total_results(i), row_units{i}, ...
        policy_E4(i), row_units{i});
end
fprintf(fileID, '\\hline\n');
fprintf(fileID, '\\end{tabular}\n');
fprintf(fileID, '\\end{table}\n');
fclose(fileID);

toc

